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A variational ground state of the repulsive Hubbard model on a square lattice is investigated 
numerically for an intermediate coupling strength (U = 8t) and for moderate sizes (from 6 x 6 to 
10 x 10). Our ansatz is clearly superior to other widely used variational wave functions. The results 
for order parameters and correlation functions provide new insight for the antiferromagnetic state 
at half filling as well as strong evidence for a superconducting phase away from half filling. 

PACS numbers: 71.10.Fd,74.20.Mn,74.72.-h 



The Hubbard model plays a key role in the analysis 
of correlated electron systems, and it is widely used for 
describing quantum antiferromagnetism, the Mott metal- 
insulator transition and, ever since Anderson's suggestion 
superconductivity in the layered cuprates. Several 
approximate techniques have been developed to deter- 
mine the various phases of the two-dimensional Hubbard 
model. For very weak coupling, the perturbative Renor- 
malization Group extracts the dominant instabilities in 
an unbiased way, namely antiferromagnetism at half fill— 
ingand <i-wave superconductivity for moderate doping 
@J3 ■ Quantum Monte Carlo simulations have been suc- 
cessful in extracting the antifcrromagnetic correlations 
at half filling 0, 0| , but in the presence of holes the nu- 
merical procedure is plagued by the fermionic minus sign 
problem [||. This problem appears to be less severe in 
dynamic cluster Monte Carlo simulations, which exhibit 
a clear tendency towards d-wave superconductivity for 
intermediate values of U 0. 

Variational techniques address directly the ground 
state and thus offer an alternative to quantum Monte 
Carlo simulations, which are limited to relatively high 
temperatures. Previous variational wave functions in- 
clude mean-field trial states from which configurations 
with doubly occupied sites are either completely elimi- 
nated (full Gutzwillerprojection) [H, 0, [l(| or at least 
partially suppressed [ll| . Recently, more sophisticated 
wave functions have been proposed, which include, be- 
sides the Gutzwiller projector, non-local operators re- 
lated to charge and spin densities 12, 0|. Our own 



variational wave function is based on the idea that for 
intermediate values of U the best ground state is a com- 
promise between the conflicting requirements of low po- 
tential energy (small double occupancy) and low kinetic 
energy (derealization). It is known that the addition of 
an operator involving the kinetic energy yields an order 
of magnitude improvement of the ground state energy 
with respect to a wave function with a Gutzwiller pro- 
jector alone ljj. In this letter, we show that such an 
additional term allows us to draw an appealing picture 
of the ground state, both at half filling and as a function 
of doping (some preliminary results have been published 



mm)- 

In its most simple form, the 2D Hubbard model is com- 
posed of two terms, H = tT + UD , with 

T = - Yl ( c L c j<?+ c ]* c ™) and L> = ^ n A n lV . (1) 

Here cJ CT creates an electron at site i with spin er, the 
summation is restricted to nearest-neighbor sites and 
rii a = c\ a Ci a . We consider a square lattice with periodic- 
antiperiodic boundary conditions and choose U to be 
equal to the bandwidth, U — St. Our ansatz 



I*) = e- hT e- 9D \<fn 



(2) 



is linked to a mean-field ground state \^o) with either a 
(d-wave) superconducting or an antifcrromagnetic order 
parameter. The operator e~ 9D partially suppresses dou- 
ble occupancy for g > 0, while e~ hT promotes both hole 
motion and kinetic exchange (close to half filling). In the 



limit h — > we recover the Gutzwiller ansatz For 
g — > oo and /i< 1 our variational problem is equivalent 
to that of the t-J Hamiltonian with respect to a fully 
Gutzwiller-projectcd mean-field state. 

The calculations for h > are carried out in momen- 
tum space where the operator D is not diagonal. There- 
fore a discrete Hubbard-Stratonovich transformation is 
applied to decouple the terms n^n^ in the operator e~ gD 
by introducing an Ising spin at each site. Expectation 
values are obtained using a Monte Carlo simulation with 
respect to the Ising spin configurations. 

We discuss first the case of an average site occupa- 
tion n = 1 (half-filling) , where the ground state is ex- 
pected to exhibit long-range antiferromagnetic order. A 
commensurate spin-density wave characterized by a gap 
parameter A^p is therefore the natural mean-field refer- 
ence state, |\&o) = |SDW). The order parameter is the 
staggered magnetization defined by 
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where N is the number of sites. In Table [T] the results ob- 
tained by minimizing the energy expectation value with 



2 



respect to the three variational parameters g, h, Aaf 
for an 8 x 8 lattice are compared with the unrestricted 
Hartree-Fock approximation (g = h = 0), the Gutzwillcr 
wave function (g > 0,h = 0) [l7j . a quantum Monte 
Carlo simulation Q and a Projector Operator technique 
[l8| . The gap parameter A^i? is very large for g = h = 
and decreases dramatically if <? and /i are optimized. We 
note that the gap parameter cannot be identified with an 
excitation gap, which in fact should vanish if a continu- 
ous symmetry is broken. The ground state energy is seen 
to vary appreciably as the parameters g and h are turned 
on and to be comparable to that found with other tech- 
niques. On the other hand, the order parameter is still 
rather large, at least in comparison to the accepted value 
of M = 0.614(1) for the 2D Heisenberg model (U = oo) 
[l9| . an upper bound for the Hubbard model. 



M 



E/t 



VMC 3.6(1) 0.89(1) -0.466(1) 

VMC 0.69 1.3 0.86(1) -0.493(3) 

VMC 3.1(1) 0.101(3) 0.32(2) 0.77(1) -0.514(1) 

QMC - - - 0.42(1) -0.48(5) 

PQ -0.521(1) 

TABLE I: Variational results (VMC) for the 2D Hubbard 
model at half-filling (8x8 lattice, U = 8t), compared to 
quantum Monte Carlo simulations (QMC) and a Projector 
Operator approach (PO). The VMC data include the unre- 
stricted Hartree-Fock approximation, the Gutzwiller ansatz 
and the present study. 

In order to extract some information about supercon- 
ducting correlations in the presence of antiferromagnetic 
long-range order, we have calculated the correlation func- 
tion Fij = (CjCj), where 
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The sites ji are the four nearest neighbors of site i and 
ctj-. = +1(— 1) in x-(y-) direction. Thus C\ creates a 
singlet pair with c?-wave symmetry centered at site i. 
is found to decrease rapidly with increasing distance, as 
expected for a gapped system. For on-site correlations 
we find F vl = 0.0637(1) for h ^ and F vl = 0.0592(1) for 
h = 0, while the results for nearest-neighbor correlations 
are = 0.0171(1) for h ^ and F 4j = 0.0155(1) for 
h = 0. The superconducting correlations are therefore 
slightly enhanced by the parameter h. 

We now discuss the effects of hole doping, and in par- 
ticular the possibility of d-wave superconductivity, as 
suggested by Renormalization Group ar gum ents Eg , 
previous variational calculations [1, lo ljiol . 11,12, 131 ] and 
quantum Monte Carlo simulations [7[. Our mean-field 
reference state is the BCS wave function with d-wave 
symmetry, \^o) = |<£BCS), characterized by a gap pa- 
rameter A describing pairing and a "chemical potential" 
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0.8125 


-0.4418(1) 


3.0(1) 


0.099(2) 


-0.849(1) 


0.8400 


-0.3972(3) 
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0.9000 


0.357(1) 
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0.106(2) 
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0.9375 
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0.110(2) 


-0.627(1) 


0.9600 


0.692(1) 


4.1(2) 


0.115(3) 


-0.583(1) 


0.9700 


0.743(1) 


4-3(2) 


0.116(3) 


-0.564(1) 



TABLE II: "Chemical potential", parameters g and h and 
total energy per site for different densities on an 8x8 lattice. 



fj, fixing the average electron density n [l5j. To reduce 
the statistical error in the Monte Carlo simulations and 
consequently the computational time, a fixed set of "Ising 
spin" configurations is first generated and then used to 
optimize the variational parameters 2^, 21 1 . 

The ground state energy and the parameters g, h, fi are 
given in Table [II] for an 8x8 lattice and various densities. 
The chemical potential \i varies strongly with doping and 
increases so much for n — > 1 that the optimization be- 
comes very difficult. The Gutzwiller parameter g also 
increases rather strongly for n — > 1, which indicates that 
the system is "more localized" at half-filling than away 
from half-filling [l(|. In contrast, the kinetic parameter 
h does not vary appreciably. 

The gap parameter A and the order parameter <f> = 



|(c-|cj.|)| are shown in Fig.QJa) as functions of the hole 
density x = 1 — n, for an 8 x 8 lattice. Both quanti- 
ties have a maximum around x — 0.1 and tend to zero 
around x = 0.18. The limiting behavior for x — > has 
not been established firmly, due to computational prob- 
lems mentioned above, but our results are consistent with 
A — > 0, $ — > 0. For the Gutzwiller wave function, the 
order parameter also exhibits a dome shape, but not so 
the gap parameter: A is found to increase monotonically 
for x — ► 0, both for finite U (inset of Fig. [2(a)) and for 
U -f oo 0. 

The condensation energy, E con d = E(0) — E(A) where 
A is the optimal gap parameter, is depicted in Fig.[T](b). 
It vanishes for x > 0.18 and increases monotonically with 
decreasing x, even beyond the hole concentration where 
both A and <I> pass through a maximum. The limiting 
behaviour for x — > is again unknown, but for x = 
antifcrromagnctism prevails. The comparison with the 
Gutzwiller wave function (inset) indicates that the addi- 
tion of the parameter h strongly enhances the conden- 
sation energy [22j. It is worthwhile to add that accord- 



ing to calculations for small clusters [14[ the difference 
AE = E var — Eq between the variational energy E var and 
the exact ground state energy E is of the same order for 
h > as the condensation energy E con d (AE « 0.007i, 
Econd ~ 0.005£ at n = 0.9), in contrast to the case h = 
where AE > E cond {AE w 0.08t, E cond w O.OOli). 

An important question is to what extent an 8 x 8 lattice 
is able to mimic the thermodynamic limit. Therefore we 
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FIG. 1: (Color online) (a): Gap (triangles) and order param- 
eter (squares) as functions of the doping for an 8x8 lattice. 
The mark at half- filling is the antiferromagnetic gap. The 
inset shows the gap parameter for the Gutzwiller wave func- 
tion, (b): Condensation energy per site. (a) + (b): Error bars 
indicate statistical uncertainties. 
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FIG. 2: (Color online) Finite-size scaling of the gap parameter 
for densities n = 0.84 (circles), n = 0.90 (triangles) and n — 
0.9375 (squares). 



have also studied other lattice sizes. The results for the 
gap parameter (Fig. [2|) show that the size effects are more 
important in regions where the gap is small (n = 0.84) 
than well inside the superconducting dome (n = 0.90 or 
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FIG. 3: (Color online) Total (circles), kinetic (squares) and 
potential (triangles) energies per site as functions of the gap 
parameter on an 8x8 lattice, for the density n = 0.9375. For 
each curve E(A = 0) has been subtracted. The relative error 
is smaller than the symbol size. The corresponding results for 
the Gutzwiller wave function are given in the inset. 



n = 0.9375). However, even at n = 0.90 where the gap 
is maximal, a 6 x 6 lattice is not large enough to give a 
reliable estimate for the thermodynamic limit. 

In Fig. [3] the kinetic and the potential energies are 
plotted separately for a density n = 0.9375 as functions 
of the gap parameter. It turns out that the maximum 
energy gain (the condensation energy) at A = 0.1 It is 
to a large extent (> 75%) due to a decrease in the ki- 
netic energy, in contrast to the BCS behaviour where 
the condensation energy is entirely due to the potential 
energy. Our findings are also qualitatively different from 
those obtained with a Gutzwiller wave function for which 
the kinetic energy increases monotonically with the gap 
parameter (inset). We have also studied the magnetic 
structure factor 

1 



i,3 



within the superconducting phase. Fig.[4]shows this func- 
tion for several densities along three different lines in the 
Brillouin zone. The structure factor is peaked at (tv,7t), 
indicating antiferromagnetic correlations. The peak de- 
creases with increasing hole concentration. The compar- 
ison with results for h = (inset) shows that the anti- 
ferromagnetic correlations are strongly enhanced by the 
parameter h. 

In summary, we have found that the addition of a 
"kinetic projector" e~ hT to the Gutzwiller wave func- 
tion yields both quantitative improvements (for instance 
for the ground state energy or for the antiferromagnetic 
order parameter) and qualitative changes (such as the 
doping dependence of the superconducting gap or the 
decrease of the kinetic energy as a function of the gap 
parameter). Nevertheless, there remains room for im- 
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In conclusion, the present variational calculations give 
an appealing picture of the ground state of the 2D Hub- 
bard model, both at half-filling and for the doped sys- 
tem. Superconductivity out of purely repulsive inter- 
actions appears very naturally in this scheme. Several 
predictions agree surprisingly well with experiments on 
layered cuprates. 
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FIG. 4: (Color online) Magnetic structure factor as a func- 
tion of the wave vector for different densities and an 8x8 lat- 
tice. The inset shows the magnetic structure factor for the 
Gutzwiller wave function at n — 0.9375. 



provement because our variational ansatz (as well as all 
trial states used previously) is linked to a delocalizcd 
mean-field reference state and thus requires a strong sup- 
pression of double occupancy, at least for U = 8t. This 
reflects the fact that this parameter regime corresponds 
to that of a doped Mott insulator which would be treated 
more naturally starting from a localized reference state. 

Finally, we comment on the relevance of our find- 
ings for layered cuprates. The antiferromagnetic ground 
state for x = and a superconducting phase with d- 
wave symmetry are well established experimentally, with 
a slightly different doping range for superconductivity 
(0.05 < x < 0.3 against < x < 0.18 in our study). The 
addition of a hopping term between next-nearest neigh- 
bor sites (parameter t') might improve this comparison. 
The typical size of the superconducting gap at optimal 
doping is 30 meV [23| , in good agreement with our result 
(A « Q.1M « 40 meV for t = 300 mcV). The com- 
parison with condensation energies extracted from spe- 
cific heat measurements is less encouraging (0.10 — 0.20 
meV/copper 0, HH against 0.005* = 1.5 meV in this 
work). Much discussion has been raised by the ques- 
tion of "kinetic energy driven superconductivity" where, 
in contrast to BCS, the energy gain arises from a de- 
crease in kinetic energy. The reported gain of kinetic 
energy AEkin ~ 0.5 — 1.0 mcV on the basis of opti- 
cal spectroscopy at optimal doping [25|, [2(| corresponds 
well to our result (AEkin ~ 1-1 meV), but we have to 
be aware that the use of a low-energy cut-off in the 
frequency-integration of the conductivity is not unam- 
biguous. Good agreement with experiment is also found 
for the strong antiferromagnetic correlations in the su- 
perconducting phase. In fact, the structure factor S(q) 
determined by neutron scattering experiments shows a 
pronounced peak at (n, n), which decreases upon doping 



[1] P. W. Anderson, Science 235, 1196 (1987). 
[2] D. Zanchi and H. J. Schutz, Phys. Rev. B 61, 13609 
(2000). 

[3] C. J. Halboth and W. Metzner, Phys. Rev. B 61, 7364 
(2000). 

[4] J. E. Hirsch, Phys. Rev. B 31, 4403 (1985). 

[5] S. R. White et al, Phys. Rev. B 40, 506 (1989). 

[6] E. Dagotto, Rev. Mod. Phys. 66, 763 (1994). 

[7] T. A. Maier, M. Jarrell, and D. J. Scalapino, Phys. Rev. 

B 74, 094513 (2006). 
[8] C. Gros, Phys. Rev. B 38, 931 (1988). 
[9] H. Yokoyama and H. Shiba, J. Phys. Soc. Jpn. 57, 2482 

(1988). 

[10] A. Paramekanti, M. Randeria, and N. Trivedi, Phys. Rev. 

B 70, 054504 (2004). 
[11] T. Giamarchi and C. Lhuillier, Phys. Rev. B 43, 12943 

(1991). 

[12] S. Sorella et al, Phys. Rev. Lett. 88, 117002 (2002). 
[13] H. Yokoyama et al, J. Phys. Soc. Jpn. 73, 1119 (2004). 
[14] H. Otsuka, J. Phys. Soc. Jpn. 61, 1645 (1991). 
[15] D. Eichenberger and D. Baeriswyl, Physica C (2007), 

doi:10.1016/j.physc. 2007.03. 366. 
[16] D. Baeriswyl, D. Eichenberger, and B. Gut, phys. stat. 

sol. (b) 244, 2299 (2007). 
[17] H. Yokoyama and H. Shiba, J. Phys. Soc. Jpn. 56, 3582 

(1987). 

[18] G. Polatsek and K. W. Becker, Phys. Rev. B 54, 1637 
(1996). 

[19] A. W. Sandvik, Phys. Rev. B 56, 11678 (1997). 

[20] D. Ceperley, G. V. Chester, and K. H. Kalos, Phys. Rev. 
B 16, 3081 (1977). 

[21] C. J. Umrigar, K. G. Wilson, and J. W. Wilkins, Phys. 
Rev. Lett. 60, 1719 (1988). 

[22] It should be mentioned that the calculations for h > 
have been carried out for a "grand canonical ensemble" 
(where jj, is fixed by the average density), whereas the 
results for h = have been obtained for a "canonical 
ensemble" (where the particle number is fixed and fi is a 
variational parameter) . The two ensembles are equivalent 
in the thermodynamic limit. 

[23] A. Damascelli, Z. Hussain, and Z.-X. Shen, Rev. Mod. 
Phys. 75, 473 (2003). 

[24] J. W. Lorain et al, Physica C 341-348, 831 (2000). 

[25] E. van Heumen et al, Phys. Rev. B 75, 054522 (2007). 

[26] H. J. A. Molegraaf et al, Science 295, 2239 (2002). 

[27] H. F. Fong et al, Phys. Rev. B 61, 14773 (2000). 



